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We propose to stabilise HMC simulations of lattice QCD with very light Wilson quarks by split- 
ting the quark determinant into two factors and by treating the factor that includes the contribution 
of the low modes of the Dirac operator as a reweighting factor In general, determinant reweight- 
ing becomes inefficient on large lattices, because the statistical fluctuations of quark determinants 
increase exponentially with the lattice volume. Random matrix theory and some numerical studies 
now suggest that the low-mode contribution to the determinant behaves differently, which allows 
factorisations to be devised that preserve the efficiency of the simulation on large lattices. 
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1. Introduction 

Simulations of lattice QCD with very light Wilson quarks are potentially affected by algorith- 
mic instabilities, sampling inefficiencies and ergodicity violations. The issue was studied in some 
detail in ref. ^ and it was concluded that simulations based on the Hybrid Monte Carlo (HMC) 
algorithm iQ] can be expected to be stable in a range of the lattice parameters and the light-quark 
masses which includes the large-volume regime of QCD at lattice spacings a < 0. 1 fm. 

The algorithm proposed here avoids the instabilities from the beginning by separating the low 
modes of the Dirac operator from the rest of the modes and by including only the latter in the HMC 
algorithm. The low modes are then taken into account by reweighting the generated representative 
ensemble of fields by the appropriate factor. Similar mode separations were recently considered by 
Jansen et al. and by Hasenfratz et al. [Q], partly with the same motivations and partly for other 
reasons (see refs. [||-^ for related earUer work). 

Low-mode reweighting tends to reduce the statistical fluctuations of observables that are sen- 
sitive to the low modes but will only work out if the reweighting factor itself does not fluctuate 
too much. From this point of view, reweighting by ratios of quark determinants (such as the ones 
considered below) does not seem to be particularly promising, because quark determinants scale 
exponentially with the volume of the lattice. Our aim in this report is to show that the situation is ac- 
tually more favourable than suspected, the main reason being that the low eigenvalues of the Dirac 
operator (except perhaps for the few lowest ones) fluctuate by no more than a distance inversely 
proportional to the lattice volume about their mean values. 

2. Determinant factorisation 

We consider lattice QCD with a doublet of light Wilson quarks and any number of heavier 
quarks. The (massive) light-quark Wilson-Dirac operator is denoted by D and the associated bare 
current-quark mass by m. On average the spectral gap of the hermitian Dirac operator 75D around 
the origin is then approximately equal to Zaw?, where Za is the renormalization constant of the 
isovector axial current [|l|]. 

Determinant reweighting starts from an exact factorisation 

det(D'''D) = W det{D'^D) (2. 1) 

of the light-quark determinant, where W is the reweighting factor and D a modified Wilson-Dirac 
operator whose determinant is included in the HMC algorithm. In the following, two choices D/, 
/ = 1 , 2, of the modified operator will be considered for which the associated reweighting factors 
are of the form 

W/ = det{w/(D"^D)} (2.2) 

(see Table 1). Note that the (complex) spectrum of 75D/ is in both cases rigorously separated from 
the origin by a distance of order /i, for all quark masses m, very much as in twisted-mass QCD with 
twisted mass /i |j8|]. The modified quark determinant det(D/Di) in fact coincides with the quark 
determinant in twisted-mass QCD. 

The inclusion of the modified quark determinant in the HMC algorithm is not expected to 
give rise to instabilities since the modified Wilson-Dirac operator is safe from having near-zero 
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Di w/(v^) w/(v^) 



2 2 

1 ^ + l-^+Olv-'*) 

Table 1: Modified Dirac operators and associated reweighting factors considered in this report. The mass 
parameter ji > can in principle be set to any value, but good reweighting efficiencies are only achieved if 
jj. is not much larger than Za»!. 

modes. Moreover, we do not foresee any difficulties in applying acceleration techniques such as 
the Schwarz preconditioning and local deflation [ [lO| | to the modified algorithm. In this report, 
however, the focus will be on the reweighting efficiency and its dependence on the lattice size. 

3. Statistical fluctuations of Wi 

We now need to distinguish the true QCD expectation value (^) of any observable ^ from its 
expectation value (^)ni in the theory with the modified quark determinant. Only the latter can be 
estimated directly using the representative ensembles of fields generated by the HMC algorithm, 
while the first is obtained through 

(3.1) 



Evidently, for the reweighting (3.1) to work out in practice, the statistical fluctuations of Wj must 
be fairly small. 

Whether this condition can be met on large lattices is unclear since 



dsi..As,Tv^{D^D + si + ... + si)-^Y (3.2) 



is the exponential of an extensive quantity X/. In particular, using the moment-cumulant transfor- 
mation one can show that 



^ = W)(Hr')=exp|£^(X-U|, (3.3) 

where {X^")con denotes the connected part of {X^"). The fluctuations of the reweighting factor thus 
grow exponentially with the lattice volume V and may therefore rapidly become too large when V 
is increased. However, as we shall see in the following sections, there are important mechanisms 
that suppress the fluctuations to the extent that determinant reweighting becomes a viable method 
in a useful range of parameters. 

4. Suppression of the high modes 

Determinant reweighting is intended for use at quark masses close to or below the range of sta- 
bility of the HMC algorithm. The values of the renormalized light-quark mass where reweighting 
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will be applied are therefore expected to be smaller than 20 MeV or so . At these quark masses, the 
low end of the spectrum of (D^D)^/^ starts at about Z^m and has an approximately constant density 
from there up to eigenvalues v„,r of at least 100 MeV. In the following, the associated eigenmodes 
of D^D will be referred to as the "low modes" of the Dirac operator and all other eigenmodes as the 
"high modes". As explained in the next section, the mass parameter /i will, in the cases of interest, 
be less than, say, 2Zpjn and therefore always well below the high eigenvalues of the Dirac operator. 

The reweighting factors W/ have been chosen so that the eigenvalues of D^D larger than 
make a monotonically decreasing contribution of order /i^'/v^' to Xi (see Table 1). Power count- 
ing then shows that the expectation values {Xf")coB are ultraviolet convergent except for {Xf)con 
which diverges logarithmically. The high modes of the Dirac operator thus contribute a term pro- 
portional to jj.'^'V to the fluctuations of Wi, with a proportionality constant that diverges at most 
logarithmically in the continuum limit. 

On current lattices the product n^V is usually much smaller than 1. For lattices of size 2L^ with 
L < 4 fm, for example, and if /Ir < 20 MeV, one obtains jX^V < 0.054. The high-mode contribution 
to the statistical fluctuations of the reweighting factors Wi thus tends to be strongly suppressed, 
particularly so in the case of W2, where there is a second suppression factor proportional to /i^. 



5. Fluctuations of the low eigenvalues 

It is still not excluded, however, that the reweighting factors Wi receive wildly fluctuating con- 
tributions from the low modes of the Dirac operator. There is no obvious suppression mechanism 
in this case, and since the number of low modes grows proportionally to V, it seems likely that 
determinant reweighting will, in practice, be limited to small lattices. 

In order to get some insight into the problem, we worked out the leading term {X^)con of 



the cumulant expansion ( |3.3[ ) in the standard two-flavour chiral random matrix theory []llp. In 
this theory, {X^)con can be expressed through the spectral density of the Dirac operator and the 



spectral 2-point correlation function, both of which are known analytically []11|,|12[]. An integration 
over the spectral parameters is then still required, but since the integrands are non-singular, it is 
straightforward to evaluate the integrals numerically. 

In random matrix theory, {X^)con is a well-defined function of the dimensionless combinations 
mLV and /iZV, where Z denotes the quark condensate in the chiral limit. To a very good approxi- 
mation, we however found that {Xf)coa only depends on the ratio of these parameters (see Figure 1). 
Random matrix theory thus suggests that the contribution of the low modes to the fluctuations of 
the reweighting factors does not change significantly with the volume V. Moreover, as can be seen 
from Figure 1 , the contribution is actually quite small up to values of n/m equal to 2 or so (note 
that Za = 1 in random matrix theory and j^/m therefore corresponds to jX^/m^ in lattice QCD). 

The outcome of our calculations in random matrix theory can be explained by noting that 
random matrices have a fairly rigid spectrum, i.e. with high probability, the low eigenvalues are 
practically unchanged from one random matrix to another. The rigidity of the spectrum may well 



'The renormalized masses are otr = Zfiju/Z^ and /Ir = ll/Z^, where Zp denotes the renormalization constant of 
the isovector pseudo-scalar density. The median v„ of the distribution of the n'th eigenvalue v„ of (Z)^£))'/^ is similarly 
renormalized through v„_r = V„/Zp. Values of these quantities quoted in physical units refer to the MS renormalization 
scheme at 2 GeV. 
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Figure 1: Values of (X,^)con computed in random matrix theory at m = 5, . . . , 30 MeV and n = 20, ... ,44 
MeV, assuming E (250 MeV)^ and V = (4.5 fm)"*^ (plot on the left). The plot on the right shows the widths 
of the distributions of the first 32 (48) eigenvalues of (D^D)'/^ on a lattice of size 48 X 24^ (64 X 32^). In 
both cases the lattice spacing and the renormalized sea-quark mass are approximately equal to 0.08 fm and 
25 MeV respectively. The grey points labelled "48 x 24^ scaled" are the 48 x 24^ data scaled by the ratio 
(24/32)"*^ of the lattice volumes. 



be related to the fact that the Vandermonde determinant, which appears in the joint distribution of 
the eigenvalues, gives rise to a repulsive force between neighbouring eigenvalues. In any case, since 
the fluctuations of the low eigenvalues are of order (LV)^^ and since there are 0{V) eigenvalues, 
their contribution to the fluctuations of the reweighting factor W/ will obviously remain bounded at 
large V. 

In lattice QCD with Wilson quarks, chiral symmetry is not exactly preserved and it is therefore 
not guaranteed that the eigenvalues of the lattice Dirac operator behave in the same way. Numerical 
studies of the 0(<3)-improved two-flavour theory however suggest that the widths of the eigenvalue 
distributions scale roughly like l/V, as in random matrix theory, except for the widths measured 
close to the threshold of the spectrum (see Figure 1)^. In the Wilson theory, the contribution of 
the low modes to the fluctuations of the reweighting factor is therefore slowly increasing with the 
lattice size, an effect that is likely to become smaller as one moves closer to the continuum limit. 



6. Computation of the reweighting factors 

An exact calculation of the reweighting factors W/ is normally not possible and actually not 
required. Stochastic estimators can be used instead, or perhaps some combination of a stochastic 
estimator and a (not necessarily exact) projector to the few lowest modes of the Dirac operator. 
Here we only consider the most obvious choice, where a. set rik{x), k = l,...,N,of pseudo-fermion 



The representative ensembles of gauge-field configurations used for our numerical studies were generated by the 



authors of ref. niBh and were made available to us through the CLS community effort. 
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Figure 2: Reweighting factor W2.24 normalized by its median at otr ~ 25 MeV and jij^/m^ = 0.7, 1.0, 1.5 
(from top to bottom), calculated for a set of independent gauge-field configurations on a 48 x 24^ (plots on 
the left) and a 64 x 32^ lattice (plots on the right) with spacing a ~ 0.08 fm. 



fields with action 



N 



^r, = £ {rik,rik) 



(6.1) 



is added to the theory and the reweighting factor W/ is replaced by 



(6.2) 



k=i 



The simulation then proceeds as before and the reweighting factor W/ a? is calculated according to 



eq. (5.2), using, for each gauge field, N randomly chosen pseudo-fermion fields. This procedure is 
correct for any A'^ > 1, but it pays to set to values significantly larger than 1, because the variance 
of Wi^N decreases when N is increased (and eventually converges to the variance of W/). 

For illustration the Monte Carlo time series of VV2.24 calculated on the lattices previously con- 
sidered are plotted in Figure 2. In all these cases, little would be gained by choosing more pseudo- 
fermion fields or by separating the lowest modes of the Dirac operator (such a mode separation 
may, however, be required at smaller quark masses). 

Figure 2 also shows that the fluctuations of W224 increase with the lattice size and that they 
are quite sensitive to the value of pL^/m^. In particular, by decreasing the latter, the fluctuations are 
quickly reduced to acceptable levels on both lattices. The fluctuations of Wi 24 at pL^/m^ = 0.7 and 
1.0 are, incidentally, practically the same as those of W2,24 at pL^/m^ = 1.0 and 1.5, respectively. On 
the lattices considered and before the performance of the HMC part of the algorithm is determined, 
it is therefore not clear whether the first or the second factorisation of the quark determinant is 
preferable. 



6 



Quark determinant reweighting 



Filippo Palombi 



7. Conclusions 

In this report we showed that determinant reweighting is likely to work out in lattice QCD if a 
factorisation of the quark determinant is chosen where the high-mode contribution to the reweight- 
ing factor is sufficiently suppressed. Somewhat surprisingly to us, the statistical fluctuations of the 
reweighting factor then do not grow rapidly with the volume of the lattice, a property that can be 
traced back to the rigidity of the spectrum of the low eigenvalues of the Dirac operator. 

We are now quite confident that simulations of the Wilson theory at very small quark masses 
can be stabilised in this way, but still need to prove this by actually performing such simulations 
using the proposed determinant factorisations. 
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